! Roel Matthysen
! s0202264
! Finale versie Project Deel 1

module methods

	implicit none
	private

	integer, parameter :: doubleprecision = kind(0.d0)
	public :: taylorexp,eig,eigqr
	
	contains

	function taylorexp(A,verbose) result(RES)
		real(doubleprecision), dimension(:,:), intent(in) :: A
		logical, intent(in) :: verbose
		real(doubleprecision) :: RES(size(A,1),size(A,2)),add(size(A,1),size(A,2))
		integer :: i
		! Add en RES worden op nul geïnitialiseerd.
		add = add*0;
		RES = RES*0;
		! Add is een eenheidsmatrix
		do i = 1, size(add,1)
			add(i,i) = 1
		enddo
		! RES wordt berekend zoals in de referentie-implementatie
		do i = 1,1000
    			RES = RES+add;
    			add = matmul(add,A)/i;
    			if(maxval(abs(add))==0) then
				call printmatrix(add,'add',verbose)
        			exit
    			endif
		enddo
	end function
	
	function eig(A,verbose) result(RES)
		real(doubleprecision), dimension(:,:), intent(in) :: A
		logical, intent(in) :: verbose
		integer :: j,i, ipvt(size(A,1)),N
		real(doubleprecision) :: RES(size(A,1),size(A,1))
		complex(doubleprecision) :: w(size(A,1)),vl(size(A,1),size(A,1))
		complex(doubleprecision) :: vr(size(A,1),size(A,1))
		N = size(A,1)
		! Eigenwaarden berekenen
		call zgeeveigen(A,w,vr)
		call printvectorc(w,'eigs',verbose)
		w = exp(w)
		vl = vl*0.0;
		! Exponentiëlen op de diagonaal
		do i=1,N
			vl(i,i)=w(i)
		enddo
		vl = transpose(matmul(vr,vl))
		vr = transpose(vr)
		! Stelsel oplossen
		call ZGESV(N,N,vr,N,ipvt,vl,N,j)
		RES = dble(transpose(vl))
	end function

	function eigqr(A,verbose) result(RES)
		real(doubleprecision), dimension(:,:), intent(in) :: A
		logical, intent(in) :: verbose
		integer :: j,i,N
		real(doubleprecision) :: RES(size(A,1),size(A,1))
		complex(doubleprecision) :: tau(size(A,1)),workqr(10*size(A,1))
		complex(doubleprecision) :: w(size(A,1)),vl(size(A,1),size(A,1))
		complex(doubleprecision) :: vr(size(A,1),size(A,1))
		complex(doubleprecision) :: eigv(size(A,1),size(A,1))
		N = size(A,1)
		! Eigenwaarden berekenen
		call zgeeveigen(A,w,eigv)
		vr = eigv
		! QR-ontbinding
		call ZGEQRF(N,N,vr,N,tau,workqr*0.d0,100*size(A,1),j)
		call printmatrixc(vr,'Result of QR',verbose)
		call printvectorc(w,'eigs',verbose)
		call printmatrixc(eigv,'eigv',verbose)
		w = exp(w)
		vl = vl*0.d0;
		! Exponentiëlen op de diagonaal
		do i=1,N
			vl(i,i)=w(i)
		enddo
		! vl blijft hetzelfde
		vl = matmul(eigv,vl)
		call printmatrixc(vl,'V*exp(eigs)',verbose)
		eigv = vr
		! Q expliciet opstellen
		call ZUNGQR(N,N,N,eigv,N,tau,workqr,5*N,j)
		call printmatrixc(eigv,'Q',verbose)
		eigv = transpose(conjg(eigv))
		! RY=Q' triangulair systeem oplossen
		call ZTRTRS('U','N','N',N,N,vr,N,eigv,N,j)
		call printmatrixc(eigv,'Y',verbose)
		call printmatrixc(matmul(vl,eigv),'imagresult',verbose)
		RES = real(matmul(vl,eigv))
		call printmatrix(RES,'final result',verbose)
		! To be deleted
	end function
	
	! Deze methodes dienen om debug-informatie naar de errorstream te schrijven
	subroutine printmatrix(A, header, verbose)
		character :: header
		logical, intent(in) :: verbose
		real(doubleprecision), dimension(:,:) :: A
		integer :: i
		if(verbose) then
		write(0,*) header
		do i=1,size(A,1)
			write(0,*) A(i,:)
		end do
		endif
	end subroutine

	! Deze methodes dienen om debug-informatie naar de errorstream te schrijven
	subroutine printmatrixc(A, header, verbose)
		character :: header
		logical, intent(in) :: verbose
		complex(doubleprecision), dimension(:,:) :: A
		integer :: i
		if(verbose) then
		write(0,*) header
		do i=1,size(A,1)
			write(0,*) A(i,:)
		end do
		endif
	end subroutine

	! Deze methodes dienen om debug-informatie naar de errorstream te schrijven
	subroutine printvectorc(A, header, verbose)
		character :: header
		logical, intent(in) :: verbose
		complex(doubleprecision), dimension(:) :: A
		integer :: i	
		if(verbose) then
		write(0,*) header	
		do i=1,size(A,1)
			write(0,*) A(i)
		end do
		endif
	end subroutine
	
	! Het enige wat deze routine doet is de eigenwaarden en eigenvectoren berekenen met behulp van DGEEV en deze in complexe variabelen steken volgens de standaard van DGEEV
	subroutine zgeeveigen(A,w,vr)
		real(doubleprecision), dimension(:,:), intent(in) :: A
		complex(doubleprecision), dimension(:) :: w
		complex(doubleprecision), dimension(:,:) :: vr
		real(doubleprecision) :: dvr(size(A,1),size(A,1)), dvl(size(A,1),size(A,1)), wr(size(A,1))
		real(doubleprecision) :: wi(size(A,1)), dwork(10*size(A,1))
		integer :: i,j, N
		logical eigvtodo
		eigvtodo = .true.
		N = size(A,1)
		call DGEEV('V','V',N,(A),N,wr,wi,dvl,N,dvr,N,dwork,100*N,j)
		do i=1,N
			if(wi(i)==0) then
				w(i) = cmplx(wr(i),0.d0,kind=doubleprecision)
				do j=1,N
					vr(j,i) = cmplx(dvr(j,i),0.d0,kind=doubleprecision)
				enddo
			else
				w(i) = cmplx(wr(i),wi(i),kind=doubleprecision)
				if(eigvtodo) then
					do j=1,N
						vr(j,i) = cmplx(dvr(j,i),dvr(j,i+1),kind=doubleprecision)
						vr(j,i+1) = cmplx(dvr(j,i),-1.d0*dvr(j,i+1),kind=doubleprecision)
					enddo
					eigvtodo = .false.
				else
					eigvtodo = .true.
				endif
			endif
		enddo
	end subroutine
	
end module
